Dark soliton states of Bose-Einstein condensates in anisotropic traps 
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Dark soliton states of Bose-Einstein condensates in harmonic traps are studied both analytically and 
computationally by the direct solution of the Gross-Pitaevskii equation in three dimensions. The 
ground and self-consistent excited states are found numerically by relaxation in imaginary time. 
The energy of a stationary soliton in a harmonic trap is shown to be independent of density and 
geometry for large numbers of atoms. Large amplitude field modulation at a frequency resonant with 
the energy of a dark soliton is found to give rise to a state with multiple vortices. The Bogoliubov 
excitation spectrum of the soliton state contains complex frequencies, which disappear for sufficiently 
small numbers of atoms or large transverse confinement. The relationship between these complex 
modes and the snake instability is investigated numerically by propagation in real time. 
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I. INTRODUCTION 

Numerous experimental studies have confirmed the 
general validity of the time-dependent Gross-Pitaevskii 
(GP) equation used to calculate the ground state 
and excitations of various Bose-Einstein condensates of 
trapped alkali atoms [|-|5|. To map the spectrum of 
collective (or particle-hole) excitations, mean- field linear 
response theories based on the Bogoliubov approxima- 
tion 1^,0 or its finite temperature extensions |^,^ have 
been developed and applied to numerous experimental 
configurations. 

The collective excitations are physically distinct from 
self- consistent excited states of the trapped gas. In the 
latter case, the stationary condensate wavefunction itself 
may contain one or more nodes. Indeed, the nonlinear 
GP equation supports many well-known self-consistent 
excitations, such as vortex states p ^| -| ri | , a nd configura- 
tions with bright and dark solitons ||l5|-|l9|| for attractive 
as well as repulsive Bose gases. These contrast quite 
strongly with the collective excitations, which are ob- 
tained from the linear response of the condensate to an 
external perturbation. 

In the case of a fundamental dark (or black) soliton, 
the condensate density vanishes along a nodal surface 
and the soliton velocity is zero. Such a solution is equiv- 
alent to two condensates with a phase difference of tt 
between them, separated by a thin impenetrable barrier, 
and is an idealization of the nodal structures obtained 
recently in a two-component system |2^. Dark optical 
solitons in nonlinear dielectric fibers have been actively 
studied pT[ | since their prediction and experimental 
observation p3|-p5[ ; the recent observation of solitons in 
trapped Bose gases has provided another striking 

manifestation of nonlinear atom optics 28 1. 

The stability of stationary dark solitons (also known as 



standing waves or kinks) in trapped condensates has been 
the subject of recent investigations jl^]. These states 
are thermodynamically unstable, since their energies are 
always higher than the nodeless self-consistent ground 
state. In addition, they can be dynamically unstable; 
in more than one dimension, an extended dark soliton 
in an optical fiber will generally undergo a 'snake defor- 
mation', where transverse modulations cause the nodal 
plane to decay into vortices ||2^. The Bogoliubov exci- 
tation spectrum for a kink is known to contain modes 
with imaginary frequencies and quasiparticle amplitudes 
localized in the notch however, the origin of these 
imaginary eigenvalues and their explicit connection to a 
dynamical snake instability remain unclear. 

In the present work, the properties and stability of self- 
consistent excited states of trapped Bose condensates are 
explored further. After a brief description in Sec. || of 
the formalism and techniques employed in the numeri- 
cal calculations, the number-dependence of the energy 
of stationary dark solitons is obtained and discussed in 
In the Thomas- Fermi (TF) limit, correspond- 



Sec. Ill 



ing to large condensates, the energy difference between 
the kink and nodeless ground states is found to be in- 
dependent of the number of atoms. In order to better 
understand this result, the soliton energy is calculated 
perturbatively aroun d the TF limit using a boundary- 
layer approach in Sec IV A . It is shown that the energy of 
the soliton state in the TF limit is identical to that of the 
'anomalous mode' in the Bogoliubov spectrum. Pertur- 
bation theory in the weakly interacting limit, carried out 
in Sec. IV B, demonstrates that this result is particular to 



large condensates, however. This perturbative approach 
also yields significant insight into the criteria for the ex- 
istence of Bogoliubov excitations with complex frequen- 
cies, discussed in Sec. IV C. The relationship between 



complex modes and dynamical instability of the kink is 
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explored in Sec. 0. In Sec. we explore the possibility 
of transferring the condensate into a kink state by a field 
excitation. The results are summarized in Sec. VII. 



II. THEORETICAL BACKGROUND 

At zero temperature, the dynamics of a single- 
component condensate are governed by the three- 
dimensional (3D) time-dependent GP equation: 



' dt 
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V2 + I4rap(r) + FH(r,i)U(r,t), (1) 



where the confining harmonic potential 

^trap(r) = i(a;2+aV + /?'^') 



(2) 



is completely anisotropic in general; in recent experi- 
ments on solitons in a Bose condensate the relevant 
parameters were a = ujy/ux = 's/S and (3 = lOz/^x — 2. 
The Hartree (mean-field) potential is written 



Vn{r,t)^Anrio\^{v,t)\'' 



(3) 



choosing the condensate wavefunction ■!/'(r, t) to be nor- 
malized to unity yields the strength parameter 779 = 
aNo/dx, where a is the atomic scattering length and iVo is 
the total number of atoms. We assume a condensate com- 
posed of Na atoms, in which case a — 52aB ~ 2.75 nm 
in Bohr radii as |2^. The above three equations are 
written in reduced units, where the length scale is dx — 
y^h/ MlOx , the time scale is T = ^'Kjujx, and the energy 
is given in units of TiiOx^ where lOx is the angular trap 
frequency in the x direction and M is the atomic mass 
of Na. 

The ground and self-consistent excited states of Bose- 
Einstein condensates are obtained by direct solution of 
the GP equation in imaginary time (t = it). At each 
imaginary time step, the chemical potential /i = {H)/Nq 
(where H is the GP operator on the right side of Eq. (|^)] 
is readjusted in order to preserve the norm of the wave- 
function. Self-consistent excitations may be found nu- 
merically by relaxation of the GP equation toward equi- 
librium, subject either to special initial conditions (spa- 
tial variations of phase or amplitude [13 30|) or applied 
constraints (such as orthogonality to the ground state). 

The initial wavefunction for the imaginary time prop- 
agation is a Gaussian V'(r, 0) = /(r) exp{— i(a;^ -|-a^y^ -I- 
P'^z^)} for small numbers of atoms Nq ^ 10^. For larger 
A^o the kinetic energy contribution to the total energy 
becomes negligible, and the initial state is chosen to 
be proportional to the Thomas-Fermi (TF) expression 

rap)/47r?7o9(^TF- 



The TF chem- 



/(i")v(MTF - Vtra 

ical potential is /itf = ^(ISa/^Tyo)^^^ in units of TiiOx and 
Q{x) is unity when x is positive and zero otherwise. The 
choice of initial state has no infiuence on the final result, 



but can improve the time required for numerical conver- 
gence. All stationary states without circulation can be 
classified by their reflection symmetry in the x, y, and z 
directions. For convenience, wavefunctions that are odd 
under a reflection in one spatial direction a = x,y, z are 
labelled Pa and are referred to as 'p-wave'. Similarly, 
the 'd-wave' states dap, with a, f3 — x,y,z, have odd 
reflection symmetry in two directions a and (3. In or- 
der to obtain p-wave or d-wave states, one may choose 
/(r) oc x,y,z or /(r) cx xy,xz,yz, respectively. Since 
the nonlinear Hamiltonian (|l|) commutes with the par- 
ity operators, however, perhaps the simplest strategy is 
to block diagonalize the GP operator according to parity 
in each direction, and solve for the lowest energy state 
in each of the eight parity manifolds. For higher lying 
stationary configurations, the solution must be explicitly 
orthogonalized to lower-energy states during the imagi- 
nary time propagation. 

We propagate the 3D time-dependent GP equa- 
tion using three distinct techniques: a variable step- 
size Runge-Kutta (RK) method second-order dif- 
ferencing (SOD) P4 I, and real-space product formula 
(RSPF) ||]. In the SOD and RK methods, the time 
propagation results from one or at most a few matrix- 
vector multiplies of the Hamiltonian onto a previously 
computed vector. The RSPF employs a split-operator 
approach which partitions the action of the kinetic energy 
matrix into a succession of 2 x 2 matrix operations on a 
known vector. It should be noted that all of these time 
propagators are efficiently implemented on distributed- 
memory massively parallel computers. 

The treatment of the kinetic energy operator forms 
the main difference in the implementation of the prop- 
agation techniques. The SOD and RSPF methods dis- 
cretize the kinetic energy operator using the simplest 
three-point, central finite difference (FD) formula for the 
second derivative operator. In the RK method, the spa- 
tial wavefunction is expanded in a discrete variable repre- 
sentation (DVR) 1^,^ based on Gauss-Hermite quadra- 
ture. The DVR has the advantage, shared also by the FD 
method, that the matrix elements of all local operators 
are diagonal and equal to their value on the spatial grid. 
The DVR kinetic energy operator is dense in each di- 
mension compared with the FD approach; however, it 
also provides a much more accurate representation of the 
derivatives than does the simple FD approximation. All 
of the methods scale formally with the number of spatial 
grid points, although the prefactor is different for each. 
Grids of the order of 200'^ points are used in the SOD and 
RSPF, while the DVR approach employs approximately 
100 functions in each spatial direction. 



III. STATIONARY STATES 

Using the diffusion form of Eq. (|l|) for the condensate 
wavefunction, the GP equation relaxes to the stationary 
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state of interest. The chemical potential /i and free en- 
ergy per particle E = fj,— ^ (Vaj/No may be obtained as a 
function of the number of atoms, and the results are sum- 
marized in Tables | and ||. The values of the chemical 
potential for the ground state agree to three significant 
figures with those reported earlier The probability 
densities in the z = plane for px and dxy dark soliton 
condensates with Nq = 2^" = 1024 are shown in Fig. |^. 
The line nodes, which are clearly visible as depressions 
in the condensate density, widen near the surface as the 
condensate density decreases and the healing length di- 
verges. 



positive energy but negative norm (or vice versa), and 
has been associated with the oscillation of a dark soliton 
in the trapped condensate at frequency w/\/2 p8||27[ |. In 
nonlinear systems, however, there is no direct relation- 
ship between the energy differences among self-consistent 
states and the collective excitations from these states; for 
example, the precession frequency (or anomalous mode) 
for an isolated vortex in a cylindrical TF condensate is 
not AEy/h, but rather ^AEy/h Indeed, as shown 
below, the energy of a black soliton displaced from the 
trap center is always smaller than hco / \f2. 
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FIG. 1. The probability densities in the 2 = 
plane for (a) and (b) d^y dark soliton states with 

= 2^° = 1024 are shown as 2D contour maps, where radial 
distances are in scaled trap units d^. Trap parameters are 

= (27r)177 rad/s, a = x/2, and /3 = 2. 

If the GP equation is relaxed subject to orthogonal- 
ity constraints with previous solutions, many additional 
stationary states may be found. The first s-wave excited 
state, with /(r) — 1 and constrained to be orthogonal 
to the ground state, has an ellipsoidal nodal surface cen- 
tered about zero. The excited state with /(r) = — 1 
and constrained to be orthogonal to the ground state, 
has two nodal surfaces which intersect the x axis. 

The numerical calculations indicate that the average 
energy per particle for a black soliton is independent of 
both geometry and particle number in the TF limit; this 
is in contrast with the energy per particle of an isolated 
vortex in a cylindrical trap, for example, which varies 
as AEy - (5/2i?2)in(i?/^) where R = and 
^ ~ 1/i? are respectively the mean TF radius and the 
healing length in units of dx ||l^. For large iVo, the en- 
ergy differences (both in the chemical potential and free 
energy) between the Px and ground states converge to 
a constant value of approximately 0.7 in units of h,LOx- 
The energy differences of the py and states are sim- 
ply scaled by a and /?, respectively; i.e. the soliton en- 
ergy is Ai?s K, 0.7 in units of htOy and hiOz- Simi- 
larly, the energies of the d^p states are approximately 
A£;, « 0.7(a-f /?). 

The energy of the black soliton relative to the ground 
state is close to the energy of the 'anomalous 

mode' in the Bogoliubov spectrum of a one-dimensional 
p-wave state in the TF limit |19|. This excitation has 



IV. PERTURBATIVE ANALYSIS 
A. Thomas-Fermi limit 

The soliton energy may be calculated using a boundary 
layer correction to the TF wavefunction ||l^,^. The px 
state requires a plane of nodes at xq = 0, but in principle 
can take any value since an oscillating dark soliton 
becomes black at its classical turning point ||l^,^ . The 
condensate density will have the TF form everywhere ex- 
cept in a small region near xo, where the kinetic energy 
drives the wavefunction to zero. 

In order to define a suitable perturbation parame- 
ter, it is convenient to further rescale the static GP 
equation (|l|), where the left side becomes /iV'(r), as fol- 
lows: {x,y,z} = R{x,y/a,z/ (3}, jl = ^i/fiUx = -R^/2, 
'7o = Voctfi/R^, and t/j"^ = tp'^af}/R^. Now, the normahza- 
tion of the wavefunction takes the form 4:TTfjo = J d^f-ip^. 
For a condensate in the TF limit with a soliton normal 
to the X-axis, one may neglect the kinetic energy con- 
tributions except in the a;-direction. The rescaled GP 
equation then effectively has cylindrical symmetry: 
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h -f + tlr 



V' = 0, 



(4) 



where e 



i? * is small in the TF limit, and = 

.x^ -I- 



In the outer (slowly-varying) region, one expands 
V'out = Xo + eXi + • • •• The TF result is recovered 
to zeroth order in e: xo — ±\/ (1 — P)/2. Of course, 
this solution is inconsistent with the boundary condition 
^{x = io) — 0; implying the existence of a boundary 
layer near x ^ xq. In this region, the outer solution is 
asymptotically xo ^ i^/ (1 — — Xq)/2. 

For the inner region it is preferable to define x = 
xq+SX, where the boundary-layer thickness is 5 <C 1 and 
|X| ^ are the regions where the inner and outer solu- 
tions must match. Since the asymptotic behavior of the 
outer solution is known , the inner wave function may be 
expanded as -ipin = ±\/{l - - xl)/2 [$o -f 5<^i + ...]. 
To lowest order in 5, Eq. (Q) becomes: 

: 92 



(1 



~x^ 



) (1 - 



$0 = 0. 



(5) 
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The 'distinguished Hmit' giving a nontrivial solution 
corresponds to S = ^/e = When \X\ > 0, 

$0 ~ 1 giving a perfect asymptotic match. With 



the substitution X = Z / ^ \{1 - - xl), Eq. (|) be- 
comes the well-known equation for a dark soliton in 
$0 + $0 - $0 = 0, yielding 



the continuum 



1 d-' 

2 Hz"- 



the exact solution for the inner wavefunction ^{X) = 
tanh X^J\{1- . 

The uniform solution for the wavefunction over all 
space may be written as V'unif = V'out + "^in — V'ovcr , where 
V'ovcr is the solution in the overlap region \X\ ^ 0: 



■0unif(p, £>) = T1 ' ^ 



l-p 



n'i — ^2 



X < tanh 



± 1 



(6) 



where the i< and i> correspond to x < xq and x > xq, 
respectively. In principle, the chemical potential may 
now be found directly from the normalization condition 
47r?7o = /rf^^V'unif- practice, however, the large num- 
ber of cross terms resulting from squaring Eq. (|^) makes 
this unnecessarily complicated. Rather, the integral over 
X is split into three regions: (I) — yl — 1^ x < Xa 

[Xa < Xo), (II) Xa < X < Xb {Xb > Xq) OT Xa < X < Xb 

{Xa — > ~oo and Xb — > oo), and (III) Xb < x < — p^- 
Since the inner and outer solutions are asymptotically 
equal in the overlap regions, the result cannot depend on 
the particular choices of Xa and Xa- One readily obtains: 



— - — (1-52)3/2 
15 6i?2^ °^ ' 



which may be inverted to yield the chemical potential for 
the soliton state 



l(15a/3%)2/^ + i=(l-i^)^/^ 
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MTF + --j={l- il) 



(8) 



The correction to the chemical potential in the TF limit 
is 1/^/2 when the soliton is at the origin xq = 0, corre- 
sponding to the p-wave state, and is zero when the soliton 
is at the surface of the cloud xq — 1. Since p = dE/dN, 
the energy per particle for the soliton is 



No 
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fJ-TF 



V2 



(1 



Xq) ftWi 



The corrections for xq — agree with the numerical 
values, and are independent of both geometry and num- 
ber at this level of approximation. Physically, the dark 
soliton has a constant energy as the number of atoms in- 
creases, because its area increases as while its width 
diminishes a,s ^ 1/R^ |l^ . A similar invariant is used 
as a measure of soliton stability in optical fibers [p^ . As 



the transverse soliton confinement becomes more appre- 
ciable (i.e. by increasing the frequencies ujy and Uz), how- 
ever, the kinetic energy in this direction will grow, and 
the low-order boundary-layer result will lose its validity. 

Two views of a black soliton state defined by Eq. (|^) 
are shown in Fig. |^. In the first, the density with a soliton 
displaced from the origin, xq = 3dx, is shown along the x- 
axis, and is compared with the TF ground state solution. 
The condensate containing a notch soliton bulges slightly 
overall in order to conserve the total number of atoms; 
the radii for the TF and soliton states are dx\/2pTF and 
dx\/2fls, respectively. In the second view, the soliton 
state is shown as a density plot in the xy-plane. The 
boundary layer theory to lowest order captures the diver- 
gence of the healing length in the vicinity of the cloud sur- 
face. In actuality, a displaced soliton would most likely 
be curved in order for the nodal plane to intersect the 
cloud boundary at a surface normal; such curvature is 
found for travelling dark solitons |27| and for displaced 
vortices in rotating trapped condensates 13 1^. 




(a) (b) 

FIG. 2. The boundary-layer approximation ^ for a con- 
densate containing a soliton at xg = Sd^ is shown for 
No = 2^** = 262,144 atoms, uj^ = (27r)177 rad/s, a = ^2, 
and /3 = 2. In (a), the view is along the a;-axis; the solid and 
dashed lines correspond to the soliton and TF ground states, 
respectively. In (b), the soliton state is depicted as a density 
plot in the xy-plane. 



B. Low density limit: Energy differences and 
anomalous mode 



The boundary layer analysis indicates that the energy 
of the self-consistent px state in the TF limit is exactly 
equal to the frequency of the anomalous mode in the 
Bogoliubov spectrum. In order to determine whether this 
holds for all densities, it is useful to consider the opposite 
limit of small condensates where analytical results can be 
easily obtained. 

The perturbation expansion for the time-independent 
GP equation (|l|) requires expanding the condensate wave- 
function ^ — tpQ + Xipl and chemical potential p = 
Po + Xpi in powers of A = 'inrjo- The normalized un- 
perturbed px state is 



■00 



1/4 



(9) 
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and jiQ = i (3 + a + /3) . Making use of the readily de- 
rived expression for a purely real condensate wavefunc- 
tion 



(10) 



one immediately obtains the first order correction to the 
chemical potential 



/ 3 /a/3 \ 



(18) 



Indeed, the perturbative corrections do not even have the 
same sign. 

C. Low density limit: Complex modes 



^ f ~ ^3 /a/3 \ 



(11) 



The low-lying excitations e may be obtained using the 
Bogoliubov equations 



Vtrap 
Vtrap 



2Hj 



zu 

-ev, 



(12) 



where u = u(r) and v = v(r). These equations may be 
written in the more convenient form (Hq + XHi) ip — e'^, 
where 



Ho = 
Hi = 







2?Ag - fli 







u 

V 



(13) 



For arbitrary trap anisotropy, there are always two de- 
generate modes with energy s = e/fnox = 1 in the un- 
perturbed Px state; these are the dipole and anomalous 
modes with positive and negative unit norms, respec- 
tively: 



^-1 = 



*2 = 



a/3 

47^3 

a/3 



1/4 



1/4 



(2x2-l)e-("'+"«'+'^^')/^ (14) 



(15) 



Note that the 'ground state' is the p-wave condensate 
wavefunction given in Eq. (^. The degeneracy between 
these states is lifted by the perturbing Hamiltonian Hi, 
and the first-order corrections to the eigenvalues follow 
directly from the diagonalization of the resulting non- 
symmetric g X g matrix 



1,2, 



,9, 



(16) 



where g is the degeneracy. Direct application of Eq. ( p^ ) 
with g — 2 yields no finite-number correction for the 
dipole mode, as expected, while the energy of the anoma- 
lous mode becomes 



It is instructive to consider the special case where all 
of the trapping frequencies ar e eq ual, a = /3 = 1. In this 
geometry, according to Ref. [[l9| , an infinitesimal con- 
densate number gives rise to pure imaginary frequencies. 
Assuming a Px state, one may define = + and 
block diagonalize the Hamiltonian into states of definite 
angular momentum Lx = mh. In the m — manifold, 
in addition to the dipole ( p^ ) and anomalous (|l^) exci- 
tations, there is a third mode with e — 1: 



(19) 



Diagonalizing the resulting 3x3 matrix (|l^), one again 
obtains no correction for the dipole mode. The remaining 
degenerate modes, however, split into complex conjugate 
pairs with energies 



1 



8V27r 



3±iV7 



hojx 



(20) 



These results are compared with numerical calculations 
for a spherical trap in Fi g. H . The numerics are extremely 
close to the expression ( P0[) for small 770, but show devi- 
ations by ?7o ^ 0.1. 



1.000 



0.995 



0.990 



0.985 

0.00 0.02 0.04 0.06 0.08 



0.015 




Tl„=iV„fl/rf„ 



(17) 



Evidently, the value of the anomalous mode in the Px 
state is not generally equal to the difference (^ in the 
chemical potentials between the p-wave and s-wave states 



FIG. 3. The real (solid lines) and imaginary (dashed lines) 
part of the complex excitation frequencies for m = are 
shown as a function of 770 = Noa/do for a spherical trap. 
Light and dark lines correspond to analytical ( po| ) and nu- 
merical calculations, respectively. 
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The corresponding complex eigenvectors both have 



zero norm J 



r [ Ui 



J d^r^'-ds^'i = 0, and 



satisfy the boundary conditions u, d ^ as r ^ cx): 

^-^ = *i - V2*2 - ^ (l - iV?) ^3; (21) 

^'^ = *i-\/2*2-^(l + j\/7)*3. (22) 

The condensate couphng between the axial and radial 
modes gives rise to modes with frequencies that are com- 
plex, rather than purely imaginary as assumed in |]l9| . It 
is important to note that the existence of complex Bogoli- 
ubov excitations does not violate the general condition on 
the quasiparticle amplitudes 0]: 



= 0. 



Hi — j and Si is complex, the corresponding particle-hole 
eigenfunction must have zero norm. 

For any given number of condensate atoms A'o in a 
cylindrically symmetric trap, there is a critical anisotropy 
ujpjujx = a such that all the Bogoliubov excitations of 
the p-wave state become purely real Indeed, in the 
limit iVo — > considered here, any a > 1 is sufficient to 
ensure the disappearance of the complex modes because 
the degeneracy between the anomalous and transverse 
modes is broken by the anisotropy. 

For a < 1, complex modes can arise in many m 
states. In this regime, there exist numerous additional 
unperturbed anomalous modes with energy e„m = 1 — 
Q;(2n + m) > but negative norm (since u = 0) that 
are degenerate with eigenstates having energy in'm' = 
a(2n' + m') — 1 = and positive norm [v — W). It is 
interesting to determine if the nonlinear coupling among 
these degenerate modes gives rise to complex excitations 
even in the limit of vanishing transverse confinement, 
a ^ 0. Consider the m = manifold and a = l/g 
with even integer q ^ 00. Degeneracies occur when the 
axial quantum numbers for both the u and v are zero and 
the radial numbers riu = ^ + p and = ^ — p, respec- 
tively, where < p < |. When p — ^, the unperturbed 
energy is eg = 1, but since the terms in Eq. (|l^) involv- 
ing cx Lg° •*/(/! will be smaller than the other terms 
by a factor 1/ql, this complex mode vanishes in the limit 
q — > 00. In the opposite limit p — 0, the unperturbed en- 
ergy is £0 = 0, and the quasiparticle amplitudes become 



9/2 



^3/4^(q/2)! ' 

Jo(^/2p)e--'/^ 



q 



(24) 



Note that this u = v solution has even x-symmetry and 
therefore does not correspond to the Goldstone mode. 
The condensate wavefunction ip in Eq. (0) becomes in- 
dependent of p, and the 2x2 matrix (Um immediately 
yields the imaginary eigenvalues 



e = ±i\l ^nQadlhuJx, 



(25) 



where no — Nq/V is the condensate density in the sys- 
tem volume V. The same solution is found for all values 
of m in this limit. Thus, even in the absence of trans- 
verse confinement, the excitation spectrum of the Px state 
contains complex modes, in agreement with the results 
ofRef. 

Additional insight into the limit of transverse decon- 
finemcnt may be gained by assuming translational in- 
variance in x and y at the outset. With the axial 
quantum numbers zero, the unperturbed energies for 
the u and v become Eq = h (/c"^ + /c"^) — 1 and Eq — 



1 — i (/cj^ + kl When these are degenerate, the off- 
diagonal couplings are non-zero only if ky=ky= ky and 
(23) k"^ = k1 ^ kz, enforcing the condition Eq = £q = and 
ky + k'^ = k"^ — 2. Note that the unperturbed energy 



£0 = and relevant wavevector |fc| = \/2 are the same as 
in the infinitely weak trap limit considered above. The 
perturbing matrix (|l^) consists of an infinite number of 
identical 2x2 submatrices for a given value of k^ and 
ky = yj2 — k^.. Each submatrix corresponds to a differ- 
ent value of m in the cylindrical case considered above, 
and yields the imaginary modes 

in agreement with the result (p5|). 

Summarizing the results of this section, in the limit of 
weak particle interactions or low condensate densities, we 
have found that the anomalous mode frequency does not 
correspond to the frequency difference between the ex- 
cited and ground state energies. In addition, solitons are 
unstable in a sufficiently loose trap. For solitons oriented 
in the radial direction of a cylindrically symmetric trap, 
imaginary eigenvalues appear in the Bogoliubov excita- 
tion spectrum when the axial frequency begins to exceed 
the radial frequency. These modes persist in the limit of 
vanishing radial confinement. 



V. COMPLEX EXCITATIONS AND SNAKE 
INSTABILITY 

The explicit connection between the existence of exci- 
tations with complex frequencies and the dynamic insta- 
bility of the p-wave state remains unclear. This 'snake 
instability is well known in the nonlinear optical com- 
munity [ pl[ , and is associated with the undulation of the 
nodal plane in the radial direction. It has been conjec- 
tured that the complex modes are responsible for the 
snake instability |19|; however, these modes have zero 



norm by definition (23), so it is not clear how they can 
become occupied. It is important, therefore, to deter- 
mine if the wavelength of the undulation matches the 
spatial dependence of complex modes in the Bogoliubov 
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spectrum, and if the snake instability disappears when 
the excitation frequencies become purely real. 

The low-lying complex modes for a stationary p-wave 
condensate are determined numerically using the Bogoli- 
ubov equations (]T^). For ease of computation, the trap 
is assumed to be cylindrically symmetric, = -\- , 
with the condensate wavefunction odd under reflection in 
axial direction x. In Fig. ^ are shown the complex modes 
in the excitation spectrum as a function of radial confine- 
ment, for a trap with axial frequency ujx = 2tt ■ 50 rad/s 
containing a condensate with 10''^ atoms. All complex 
modes have even axial parity. For relatively weak trap 
anisotropy lOpjiOx = a ^ 3, there are several complex 
modes with angular momentum projections to = 0, 1, 2. 
As a increases, these modes disappear in turn until only 
one purely imaginary mode with m = 1 remains. Its 
magnitude reaches a maximum at a w 6, and vanishes at 
a « 10. A similar removal of the imaginary modes can 
be effected by decreasing the particle number at fixed 
geometry. 

2.0 I — ' — ^ — ' — ^ — ' — ^ — ' — ^ — ' — ^ — ' — ^ — ' — ^ — ' — ^ — ' — I 




FIG. 4. The real (filled symbols) and imaginary (open sym- 
bols) parts of the low-lying complex excitation frequencies are 
given as a function of the trap anisotropy a = cjp/ujx for a 
cylindrically symmetric trap, where the condensate is in a px 
state with 10'* atoms and uj^ ~ 27r-50 rad/s. Excitations with 
m = 0, 1, and 2 are represented by circles, squares, and dia- 
monds, respectively; the dashed line denotes the TF estimate 
of the anomalous mode. 

For large anisotropy a 3> 1, the stability criterion is 
expected to be approximately Oc > m/2.4 [|19|, where 
fl — ^/hujx- In the TF approximation, ac > 7.6 for the 
geometry considered above. The larger value required 
here is due to deviations from the TF limit; the radial 
wavefunction approaches a Gaussian when the transverse 
confinement is strong. When a = 10, the TF chemical 
potential is /ixp — 23.88?iWa: while the actual value is 



determined to be « 2Q.'&^hwx- 

The dynamic stability of dark soliton excited states is 
investigated by propagation of the GP equation in real 
time for an extended period. In the absence of an applied 
perturbation, the self-consistent states should remain ab- 
solutely stationary. In practice, however, numerical noise 
inherent in the propagation algorithm is magnified by the 
nonlinearity. Although the norm and chemical potential 
are conserved to one part in 10^^ and 10^^ respectively 
during the numerical propagation, the kink state eventu- 
ally decays. In order to make contact with the complex 
excitation frequencies in Fig. ^, we considered the cases 
a = 1, 6, and 10. For the first two cases, the modes 
with the largest imaginary component are e — Q.Q'iifujJx 
and l.TAiTiLOx', the p-wave states were found to decay in 
real time in approximately 160 ms and 40 ms, respec- 
tively. The third case with a = 10 remained stable for 
the longest propagation time considered, 200 ms. Thus, 
the lifetime of the p-wave state appears qualitatively to 
scale with the inverse of the largest imaginary mode in 
the excitation spectrum. 
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FIG. 5. The spatial variation of the mode in the Bogoliubov 
spectrum with the largest imaginary component e = 2.16ihujx 
is shown for a px-wave condensate containing 10^ atoms in a 
spherical trap ujipjiOx = 1 with = Stt ■ 50 rad/s. The axial 
(p = 0, a; > 0) and radial (x = 0) dependences of the complex 
u amplitude are shown in (a) and (b), respectively (note that 
i]} and w are odd and even in x, respectively, and that |u| = \v\ 
for this pure imaginary mode). The axial dependence of the 
condensate wavefunction is shown for comparison in (a). 

As illustrated in Figs. || and ^, there is a close similar- 
ity between the spatial variation of the cigcnmodc with 
the largest imaginary component and that of the soliton 
nodal plane during the initial decay. For a p-wave state 
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in a sperical trap, such as the one considered above but 
with 10^ atoms, the relevant excitation is purely imagi- 
nary with an energy e = 2.16ihujx, and \u\ = \v\ (as is the 
case for all pure imaginary modes). Fig. ^ shows the cor- 
responding radial and axial dependences of the complex 
u Bogoliubov amplitude. The quasiparticle amplitudes 
are highly localized axially, but oscillate radially within 
the soliton nodal plane. It is interesting to note that the 
imaginary components of the excitation energies tend to 
decrease as the number of radial nodes increases; this 
behavior is due to the effective negative kinetic energy 
of the dark soliton and is reminiscent of internal 
waves at the interface between layers in stratified fluid 
mixtures iSSl. 
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(c) 
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all shape follows closely that of the largest u amplitude 
shown in Fig. |^: the two radial nodes of this imaginary 
excitation correspond to stationary points in the soliton 
bending. At an intermediate time, Fig. |^(e), the soliton 
has decayed into two concentric vortex rings whose cores 
are located at these nodes. The outer vortex ring decays 
rapidly (by 75 ms) to the condensate surface, where it 
shrinks considerably, as shown in Fig. |^(f); the inner ring 
was found to remain stable for much longer times. The 
vortex rings are barely visible in the integrated densities 
shown in Fig. ^a)-(c), so we expect the experimental 
observation of these intriguing features using standard 
absorption imaging to be a challenge. The decay of the 
soliton into vortex rings produces a large number of den- 
sity oscillations, which are required in order to conserve 
the total energy of the system and which are undamped 
in the present formalism. At long times, these oscilla- 
tions are most evident at the condensate surface, giving 
rise to the bright halo in Fig. ^(f). 



(d) 



(e) 



(f) 



FIG. 6. Snapshots of the snake instability are shown for a 
Py-wave condensate containing 10^ atoms in a spherical trap 
with oj = 27r ■ 50 rad/s. Times after the initial formation 
of the soliton state are 47 ms, 50 ms, and 77 ms for (a)-(c) 
and (d)-(f). In (a)-(c), the brightness is proportional to the 
condensate density, and the images correspond to densities 
integrated down the line of sight. In (d)-(f), the brightness is 
inversely proportional to the condensate density, and regions 
outside the TP sphere are rendered transparent in order to vi- 
sualize nodes in the condensate interior; the color corresponds 
to the phase: = through 2n is represented by the sequence 
red-green-blue-red. The view is perpendicular to nodal plane; 
prior to the snake instability the black soliton would appear 
as a featureless disk. 

In Fig. H, snaphots of the snake instability are shown 
for a Py-wave state in a spherical trap [for convenience, 
the axis perpendicular to the nodal line is taken to be 
along y, the vertical direction in Fig. ^(a)-(c)]. The GP 
equation is solved on a Cartesian mesh with no parity re- 
strictions. After approximately 40 ms of real time propa- 
gation, the black soliton begins to undulate. The spatial 
variations are symmetric about the y-axis, originating 
near the center and propagating outwards. The over- 



FIG. 7. The breakup of a py state is shown as a function 
of time for A'^o ~ 10® atoms, lOx ~ (27r)14 rad/s, a = y/2, 
and /3 = 2. From the top left to the bottom right in raster 
order are shown times t = 15 ms through 20 ms in 1 ms 
increments after the initial state is formed. The view is along 
y, and the Hamiltonian was constrained to even parity along 
X and z for ease of computation. The rendering is identical 
to that of Figs, ^(d)-(f). The filamentation is almost entirely 
constrained to the original nodal xz-plane. 
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It is interesting to investigate the snake instability 
for parameters relevant to recent experiments on dark 
solitons in trapped condensates |2^. Fig. |^ shows the 
breakup of a black soliton in a completely anisotropic 
trap containing one million Na atoms. The undulations 
are already pronounced by 12 ms, and radiate outwards 
from the soliton center as found above. Unlike the spher- 
ically symmetric condensate shown in Fig. however, 
the soliton does not decay into concentric vortex rings 
far from its center, but rather into a series of approxi- 
mately evenly spaced curved vortex lines in the direction 
of weakest confinement (taken to be x). The results im- 
ply that the largest imaginary mode in the Bogoliubov 
excitation spectrum has an energy of ~ 6hu)x and has 14 
nodes along x. At longer times, the innermost two vortex 
rings make contact with one another, and subsequently 
detach into a vortex line and a ring. The simulations in- 
dicate a rich dynamics among quantized vortices in these 
systems that are only beginning to be explored js^ . 

VI. FIELD EXCITATION 

Several techniques have been proposed for the exper- 
imental production of dark solitons in trapped conden- 
sates, including adiabatic Raman transitions to the p- 
wave state |l^, preparation of the condensate in a su- 
perposition of two internal states ||l^ , collisions between 
two separate condensates , and phase imprinting ; 
the last approach was recently implemented experimen- 
tally |2^,|2^. For small particle interactions, the self- 
consistent excited states (such as a p-wave or d-wave 
condensate with at least one stationary dark soliton) ap- 
proach the non-interacting single-particle excitations of 
the harmonic trap. In this regime, it might be possible 
to transfer most of the condensate into a soliton state 
by applying an external field resonant with the energy 
difference between the ground and excited self-consistent 
states. 

We have attempted to excite the dxy dark soliton state 
for a small condensate containing A^o = 1024 atoms in a 
completely anisotropic trap, shown in Fig. ^b). The 
field excitation is modeled by a large-amplitude time- 
dependent spatial perturbation in the GP equation (|^) 
given by T4dd(i",i) = A{t)xy cos{cjpt), where the ampli- 
tude A{t) is 25% of a trap energy huj^ and includes a 
smooth turn on and turn off as a function of time. The 
probe frequency is set at ujp = 2.06, which is the sepa- 
ration in chemical potential of the ground and first dxy 
state (cf. Table H). 

The time-dependent probability density in the z — 
plane is shown in Fig. g with snapshots at t — Tp, 4Tp, 
7Tp, and lOT^, where Tp = 2n/LUp w 0.5T. At t = Tp, 
the wavefunction overlap with the ground state is 0.98 
and with the dxy state is 0.01. Until t = 4Tp, the wave- 
function overlap with the ground state decreases mono- 
tonically to 0.57, while the overlap with the dxy state 



increases to a maximum of 0.25. Thereafter, as illus- 
trated at times t = 7Tp and t = lOTp, the wavefunction 
overlaps with both the ground and this particular dark 
soliton state decrease. While the perturbing potential is 
not able to yield a dxy state, the last image in Fig. ^ at 
t = lOTp shows the formation of a number of depressions 
in the probability density at the very edge of the conden- 
sate, corresponding to multiple vortices with both senses 
of circulation. 




5.00 2.50 0.00 2.S0 5.00 5.00 2.50 0.00 2.50 5.00 



FIG. 8. The probability density in the z = plane for the 
time evolution of the ground state with A'^ = 1024 under a 
time- varying spatial perturbation is shown as a contour map 
at (a) t = Tp « O.ST, (b) t = 4Tp, (c) t = 7Tp, and (d) 
t = lOTp. Radial distances are in scaled trap units dx, and 
trap parameters are ujx = (27r)177 rad/s, a = \/2, and P — 2. 
In (d), density minima at the top and bottom correspond to 
vortices, those at left and right to antivortices. 

The inability of a field excitation to produce dark soli- 
tons in not surprising, in light of the snake instabilities 
discussed in Sec. ^ Dark solitons are unstable against 
the formation of vortices under ideal conditions, and are 
increasingly unstable in the presence of an external per- 
turbation. Nevertheless, applying a very large amplitude 
perturbation may be a viable alternative approach for 
the formation of vortices in systems with much larger 
condensate densities; this possibility remains for future 
work. 



VII. CONCLUSIONS 

In summary, self-consistent ground and excited states 
of a condensate in an anisotropic harmonic trap are calcu- 
lated by direct solution of the time-dependent GP equa- 
tion in three dimensions. The energy of a dark soliton at 
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the center of the trap, relative to the true ground state, 
is found to be independent of both the number of atoms 
and the trap geometry in the TF hmit, with a value equal 
to the soliton oscillation frequency. In the weakly inter- 
acting limit, however, the energy of the anomalous mode 
does not equal the energy of the dark soliton. In both 
limits, the low-lying Bogoliubov excitation spectrum of 
p-wave states is found to contain modes with complex 
frequencies, which may be removed by strong trap con- 
finement in the soliton plane or by decreasing the con- 
densate density. These complex modes are shown to give 
rise to the snake instability of the solitons observed in 
real time propagation of the GP equation. In extended 
trap geometries, the solitons are found to decay into vor- 
tex lines and rings at long times, imposing constraints 
on the ability to generate and observe these intriguing 
excitations in current experimental geometries. 
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TABLE I. The chemical potential /i and free energy per 
particle E of the ground state (subscript 0) and p-wavc dark 
solitons oriented along x, y, and z are given as a function of 
the number of atoms in the condensate A^o in units of hcjx- 



No 


no 


Eo 


Hx 


Ex 


Hy 


Ey 




E, 


2io 


3.57 


2.99 


4.40 


3.87 


4.79 


4.27 


5.36 


4.85 


2l2 


5.42 


4.20 


6.20 


5.02 


6.56 


5.40 


7.09 


5.95 


2l4 


8.90 


6.59 


9.64 


7.36 


9.98 


7.71 


10.47 


8.23 


2l6 


15.13 


10.96 


15.85 


11.70 


16.17 


12.03 


16.63 


12.52 


2i8 


26.10 


18.75 


26.80 


19.47 


27.12 


19.79 


27.57 


20.25 


220 


45.28 


32.41 


45.94 


33.10 


46.28 


33.43 


46.72 


33.88 



TABLE II. The chemical potential ji and free energy per 
particle E of the ground state (subscript 0) and d-wave dark 

solitons with nodes along (x, y), (x, z), and (y, z) are given as 
a function of the number of atoms in the condensate No in 
units of TvjJx- 



No 


no 


Eo 


t^x,y 






Ex,z 


f^y,!: 




2io 


3.57 


2.99 


5.63 


5.16 


6.20 


5.74 


6.60 


6.15 


2l2 


5.42 


4.20 


7.33 


6.22 


7.87 


6.77 


8.23 


7.16 


2l4 


8.90 


6.59 


10.72 


8.48 


11.21 


9.00 


11.55 


9.35 


2i6 


15.13 


10.96 


16.89 


12.77 


17.35 


13.26 


17.67 


13.59 


2i8 


26.10 


18.75 


27.83 


20.51 


28.27 


20.97 


28.59 


21.29 


220 


45.28 


32.41 


46.94 


34.12 


47.38 


34.57 


47.73 


34.90 
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